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We study the four-point dynamic susceptibility X4(i) obtained from Brownian dynamics computer 
simulations of the Kob- Andersen Lennard- Jones mixture. We compare the results of the simulations 
with qualitative predictions of the mode-coupling theory. In addition, we test an estimate of the 
four-point susceptibility recently proposed by Berthier et al. [Science, 310, 1797 (2005)]. 
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I. INTRODUCTION 

The origin of the extreme slowing down of liquids' dy- 
namics upon approaching the glass transition and the 
very nature of this transition are still hotly debated de- 
spite many experimental, simulational, and theoretical 
studies performed in the last two decades One of 
the fundamental difSculties is that the slowing down ap- 
peared to be a local, small scale phenomenon which is 
not accompanied by a growing correlation length. No 
long-range correlations have ever been found in any static 
quantity upon approaching the glass transition. The 
attention has recently shifted to dynamic correlations 
0, and there is evidence that there is a dynamic cor- 
relation length that slowly grows upon approaching the 
glass transition. Unfortunately this dynamic correlation 
length cannot be easily obtained from experimental data. 

The dynamic correlation length is often obtained from 
the four-point dynamic susceptibility Xii^)- This sus- 
ceptibility is related to a space integral of a four-point 
density correlation function that quantifies correlations 
between relaxation processes at different points in space. 
It is assumed that the increase of the dynamic suscep- 
tibility signals the growth of the range of the correla- 
tions between relaxation processes. Furthermore, the 
four-point dynamic susceptibility can be used to discrim- 
inate between different theoretical approaches to glassy 
dynamics. For example, in a recent paper (3| Toninelli et 
al. discussed predictions for X4(0 obtained from various 
theoretical approaches to glassy dynamics and compared 
these predictions with two different simulations of atom- 
istic models of glass forming liquids. 

A problem with the four-point susceptibility is that it 
is difficult to extract from experiments. This problem 
has been addressed in another recent paper . Berthier 
et al. argued that the four-point susceptibility can be 
estimated using derivatives of an intermediate scattering 
function with respect to thermodynamic parameters such 
as temperature or density. Since the intermediate scat- 
teringfunction can be easily obtained from experiments, 
Ref. |j| introduced a way to experimentally investigate 
the existence of a growing correlation length. 

The goal of this contribution is twofold. First, we an- 
alyze the four-point susceptibility of the Kob- Andersen 
[5| Lennard-Jones binary mixture undergoing Brownian 



dynamics. Second, we test the approximate estimate for 
the four point susceptibility proposed by Berthier et al. 

A- 

The paper is organized as follows: in Sec. ^l^e briefly 
review the details of the simulation; in Sec. 11111 we de- 
fine the four-point susceptibility and analyze the simu- 
lation results for this quantity; in Sec. IIVI we compare 
the four-point susceptibility obtained directly from com- 
puter simulations to the approximate expression derived 
by Berthier et al.; and in Sec. 0we discuss our findings. 



II. 



SIMULATION DETAILS 



We simulated a binary mixture that was introduced by 
Kob and Andersen which consists of Na — 800 particles 
of type A and Nb = 200 particles of type B. The inter- 
action potential is Vafj{r) = ^eai3[{<yai3/TY'^ — {<7a()/r)^]i 
where a,/? e {A,B}, caa = 1-0, aAA = 1-0, cab = 1-5, 
Cab — 0.8, esB = 0.5, and <jbb = 0.88. The interaction 
potential was cut off at 2.5 aap. We used a cubic simu- 
lation cell with the box length of 9.4 aAA with periodic 
boundary conditions. 

We performed Brownian dynamics simulations. The 
equation of motion for the position of the ith particle of 
type a, rf, is 



t.0 



(1) 



where is the friction coefficient of an isolated particle, 
^0 = 1-0, and F°' is the force acting on the ith particle of 
type a, 



3 /3=1 



(2) 



with Vf being the gradient operator with respect to ff 
(note that the term with (3 = a and i — j has to be 
excluded from the double sum in Eq. (O). In Eq. 
the random noise ffi satisfies the fiuctuation-dissipation 
theorem. 



(3) 



Here Dq is the diffusion coefficient of an isolated particle, 
Dq = fes^/^o, where fcs is Boltzmann's constant. Fur- 
thermore, in Eq. (j2Jl 1 is the unit tensor. The equations 
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of motion H1I3() allow for diffusive motion of the center 
of mass, thus all the results will be presented relative to 
the center of mass {i.e. momentary positions of all the 
particles are always relative to the momentary position 
of the center of mass 0). The results are presented in 
terms of the reduced units with aAA, ^aa, ^AA/ks, and 
'^AA^o/^AA being the units of length, energy, tempera- 
ture, and time, respectively. In these units the short-time 
self-diffusion coefficient is proportional to the tempera- 
ture, thus the time is rescaled to a reduced time equal 
to tDo/a^^ to facilitate comparisons with theoretical 
approaches that often assume temperature-independent 
short-time dynamics. 

The equations of motion were solved using a Heun al- 
gorithm with a small time step of 5 x 10^^. Wc simulated 
a broad range of temperatures 0.44 < T < 5.0. Here we 
present results for the following temperatures: T — 0.45, 
0.47, 0.50, 0.55, 0.60, 0.65, 0.80, and 1.0. We ran equi- 
libration runs and 4-6 production runs. The equilibra- 
tion runs were typically twice to four times shorter than 
the production runs, and the latter were up to 1.2 x 10^ 
time steps long for the lowest temperature discussed here, 
T = 0.45. The results presented are averages over the 
production runs. 



III. FOUR-POINT DYNAMIC 
SUSCEPTIBILITY 

In the context of the dynamics of supercooled liq- 
uids, the four-point susceptibility has been introduced 
by Glotzer and collaborators 1^. The four-point suscep- 
tibility that we discuss in this paper is slightly different 
from that defined in Refs. 0- We consider here the sus- 
ceptibility which was extensively analyzed by Toninelli et 
al. 3] , which is defined as the variance of the fluctuations 
of the self-intermediate scattering function. 

We start with the definition of the self-intermediate 
scattering function, Fs{k;t), 

Fs{k;t) = ^l^cosfc. [rl(t) -rl(0)]^ . (4) 

Next, we define the fluctuation of the instantaneous value 
of the scattering function, 6Fs{k;t), 

SF.ik; t) = 1 ^ cos fc • [n{t) - niO)] - F,ik; t). (5) 

i 

The four-point susceptibility, X4(^)) then defined as 

X4it) = N (^SF,{k;t)dFs{k;t)y (6) 

The four-point susceptibility depends on the wave vector 
k (and it is sometimes denoted by Xfc(i))- This wave vec- 
tor is customarily fixed at theposition of the maximum 
of the static structure factor m. 



The system considered in this paper is a two- 
component mixture, thus instead of one self-intermediate 
scattering function Fs{k] t) we could introduce two differ- 
ent scattering functions involving particles A and B. All 
the results presented in this paper concern the A particles 
only, thus, e.g. sums in Eqs. (|4I5|I run only over the A 
particles. Since we are not presenting any results for the 
B particles we do not introduce additional sub-scripts or 
super-scripts indicating particle labels. The magnitude 
of the wave vector k is fixed at the position of the max- 
imum of the partial structure factor of the A particles, 
|fc| = 7.25. 

The mode coupling theory of the glass transition was 
formulated by Gotze and collaborators for Newtonian 
systems , and was later extended by Szamel and Lowen 
to Brownian systems 'lO^ . The theory makes predictions 
for the self- and collective intermediate scattering func- 
tions. Biroli and Bouchaud jfj] recently argued that 
the mode coupling theory could be understood as a dy- 
namic mean-field theory, and that the usual mode cou- 
pling equations are saddle point equations obtained from 
an action functional. This new interpretation made it 
possible to calculate fluctuations of the order parame- 
ter (i.e. fluctuations of a two-point dynamic correlation 
function) from the inverse of the second derivative of the 
action functional. The details of the calculation have not 
been published, but the main predictions have already 
been discussed and compared with Newtonian dynamics 
simulations 3]. According to Biroli and Bouchaud, on 
the [3 relaxation time scale the four-point susceptibility 
grows with time as a power law, X^it) oc ^ with expo- 
nents equal to the standard mode coupling exponents a 
and b in the early and late (3 regime, respectively. Fur- 
thermore, Xi{i) reaches its maximum value, X4(i*), on 
the time scale of the a relaxation time, t* ~ Tq,, and 
the maximum value diverges upon approaching the mode 
coupling temperature, X4(^*) oc (T — Tc) ''^^ with 71 = 1. 

The authors of Refs. 0,13 emphasized that these pre- 
dictions are valid for Newtonian systems in the micro- 
canonical (NVE) ensemble and hinted that they may be 
different in the canonical (NVT) ensemble. We present 
results obtained from Brownian dynamics simulations. 
In such simulations constant temperature is maintained 
automatically by the equations of motion. Thus, in prin- 
ciple it is not clear whether the predictions of Biroli and 
Bouchaud [Til are relevant for our simulation results. 
However, there are subtle theoretical arguments that sug- 
gest that the reverse is true 0. In addition, at the 
level of the mode coupling equations, there is no differ- 
ence between Newtonian systems in the NVE or NVT en- 
semble, and Brownian systems (beyond the microscopic 
time scale, i.e. on the time scales of the /3 relaxation and 
longer) 10] . These reasons encouraged us to compare our 
results obtained from Brownian dynamics simulations to 
the predictions of Biroli and Bouchaud . 

In Fig. ^we show the general shape of the four-point 
susceptibility, X4(t), for several temperatures. The time 
dependence of X4(i) is similar to that obtained for Newto- 
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FIG. 1: Time dependence of the four-point susceptibility 
X4{t) plotted vs. tDo/a'^ for T =1.0, 0.80, 0.60, 0.55, 0.50, 
0.47, and 0.45 (left to right). 



nian systems, see Figs. 4 and 5 of Ref. In particular, 
it can be argued that there is power law like growth in 
time of the four-point susceptibility as it approaches it's 
maximum value. This power law-like dependence will be 
further analyzed at the end of this section. 

In Fig. [3 we present the results of some quantitative 
analysis of Xii^)- Fhst, in Fig. |3i we show the tem- 
perature dependence of the its maximum value, Xii^*)- 
We plot the maximum value vs. e = (T — Tc)/Tc where 
Tc ~ 0.435 is the standard value of the mode coupling 
temperature for the Kob- Andersen Lennard-Jones binary 
mixture j^. In an intermediate range of temperatures, 
X4(i*) grows as a power law with decreasing {T — Tc)/Tc, 
X4(i*) OC e~'^i. The exponent obtained from the fit, 
7i — 0.995 ± 0.05, is very close to the theoretical pre- 
diction of Biroli and Bouchaud, 7^'' = 1. It should be 
noted that the range of reduced temperatures for which 
the power law dependence of X4(i*) is observed coincides 
with the range of reduced temperatures for which mode 
coupling theory was found to correctly describe the time 
evolution of the self-intermediate scattering function and 
the mean square displacement |13| . In Fig. we com- 
pare the temperature dependence of the time at which 
Xi(t) reaches its maximum value, i*, with that of the a 
relaxation time, Tq,. The latter time is defined as the time 
at which the self-intermediate scattering function decays 
to of its initial value, Tq.) — e^^. We find that i* 
and Ta are very close and have the same temperature de- 
pendence. In particular, in the same intermediate range 
of temperatures these times grow according to power laws 
with decreasing (T — Tc)/Tc, t* oc e~"' , Tq oc e~'^. The 
exponents obtained from the fits, 7* = 2.27 ± 0.04 and 
7 = 2.31 ± 0.02 are very close. The numerical solution 
of the mode coupling equations predicts a slightl y h igher 
value of the scaling exponent for r^, 7*^* — 2.46 [l3- Fi- 
nally, in Fig.Ht we plot Xiit*) vs. t* . We find that shght 
deviations from power laws that are visible in Figs. 
and |2h are magnified in Fig. Nevertheless we can 
still fit a power law Xiit*) ^ (i*)^/^^ with the exponent 




FIG. 2: (a) Temperature dependence of the maximum value 
of the four-point susceptibility: X4(^*) plotted vs. e = 
(T - Tc)/Tc. Solid line is the power law fit X4(i*) oc e"^i 
with 71=0.995. (b) Temperature dependence of the time at 
which X4(^) reaches its maximum value, t* (circles), and the 
a relaxation time, Tq (squares). Both times are plotted vs. 
e = (T — Tc) /Tc- The solid line and the dashed line are power 
law fits t* oc with 7* =2.27 and Tq oc with 7=2.31 

, respectively, (c) The maximum value of the four-point sus- 
ceptibility, X4(^*)i plotted vs. t* . The solid line is the power 
law fit Xiif) « (t*)^''^" with 72 = 2.31. 



72 = 2.31 ±0.21. 

In Fig. 131 we address the question of the power law de- 
pendence of the four-point susceptibility on time. Instead 
of trying to fit power laws to Xi^^) over some specific 
time intervals, we have numerically calculated the deriva- 
tive of lnx4(i) with respect to Int, c?lnx4(t)/rflnt = 
(t/x4{t))dx4(t)/dt. For any time interval over which 
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FIG. 3: Derivative of Inxiit) w.r.t \nt, d\nx4(t)/d^nt — 
{t/X4{t))dX4{t)/dt, plotted vs. t/t* for T = 0.45, 0.47, 0.50, 
and 0.55 (left to right). 



X4 (i ) depends on f in a power law fashion, d In ^4 (t ) /d In t 
should have a constant value (a plateau). On the basis 
of the predictions of Ref. 0| we expect two different 
plateaus in our data. A plateau in the early-/3 regime 
and another one in the late-/? regime, i.e. on approaching 
the maximum of X4(i) (additional plateaus are predicted 
for shorter times, and we do not analyze them here). 
However, the simulations agree with mode coupling pre- 
dictions only over a restricted range of temperatures. It 
is not clear whether these two time scales are well sepa- 
rated over this temperature range, thus it is not clear if 
we could see two different plateaus. 

Indeed, Fig. |31does not show well developed plateaus. 
However, we notice that Xi{t) grows with t with an ex- 
ponent of approximately 6 = 0.8 upon approaching it's 
maximum for all temperatures. This exponent is compa- 
rable to the exponents obtained from Newtonian dynam- 
ics simulations Q ■ It is quite a bit higher than the expo- 
nent obtained from the mode coupling theory 6*'' = 0.62 

M. 



IV. TEST OF THE ESTIMATE OF THE 
FOUR-POINT SUSCEPTIBILITY PROPOSED BY 
BERTHIER ET AL. 

The four-point susceptibility can be easily obtained 
from computer simulations, but it is not readily de- 
termined from experiments. To address this problem 
Berthier et al. proposed an approximate estimate for 
the four-point susceptibility in terms of the derivative of 
the self-intermediate scattering function with respect to 
temperature, XT{t), 



XT{t) 



dFsik;t) 
dT ■ 



(7) 



The starting point of their argument was a fluctuation- 
dissipation relation (Eq. (2) of Ref. P|). If we naively 
adopted this relation to our system [i.e. Brownian dy- 



namics, NVT ensemble) we would get 

kBT^XT{t) - l5Fs{k-t)6V{Q) 



(8) 



Here 5V{t) denotes the instantaneous fluctuation of the 
total potential energy. 



i,j al3 



SVit) = Vit) - {V{t)) 



(9) 



(10) 



The derivative with respect to the temperature in 
Eq. HIV|I . Eq. ||SJ) and in the remainder of this section 
have to be calculated while keeping the short-time diffu- 
sion coefficient Dq constant since our short-time diffusion 
coefficient is proportional to the temperature. 

The main analytical result obtained by Berthier et al. 
was an exact lower bound for Xiit) in terms of XT(t) 
(Eq. (5) of Ref. ,4]). Naively applying this result to our 
Brownian system we would get 



4 



(11) 



Here c^°' is the potential contribution to the constant 
volume specific heat per particle. 

Berthier et al. found that the difference between the 
right and left sides of relation (|ll|l diminishes with de- 
creasing temperature. On this basis, they proposed using 
the right-hand-side of Eq. (|ll|l as an approximate esti- 
mate for the four-point susceptibility. Furthermore, they 
showed that this estimate can be easily calculated using 
either experimental or simulational results for the self- 
intermediate scattering function. 

Both the fluctuation-dissipation relation and the 
bound derived by Berthier et al. are valid only if the 
equations of motion do not involve the temperature. The 
only place where the temperature can enter is the initial 
condition which is given by the canonical distribution. 
This is true neither for the usual NVT computer simula- 
tions, where the temperature enters into the equations of 
motion via a thermostat, nor for our Brownian dynamics 
simulations, where the temperature enters via the noise 
strength. Thus we have no arguments in favor of either 
the fluctuation-dissipation relation ((HJ or the bound 1111) 
for our system. In fact, we show explicitly in Fig. 2|that 
the relation (jSJ is violated in Brownian systems. On the 
other hand we do not have strong numerical evidence for 
the violation of the inequality lfTT|l . Figure suggests 
that this inequality is violated, but the extent of the vi- 
olation is smaller than the error bars. 

We should mention that there is a related, exact bound 
for the four-point susceptibility that follows from the 
Cauchy-Schwarz inequality (see footnote 22 of Ref. 0): 



Xiit) > 



SFsik;t)6ViO) 



k„T^d 



pot 



(12) 
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FIG. 4; Comparison of fcsT\T(t) = ksT'^ dF{k;t) /dT 
(closed circles) with (^5Fs(k;t)5V{Q)^ (open squares) at three 
representative temperatures: (a) T=1.00, (b) T=0.55 and (c) 
T=0.45. 



However, as we show in Fig. |21this bound does not lead 
to a useful estimate for the four-point susceptibility. 

Berthier et al. also gave a different, general argument 
leading to an estimate for the four-point susceptibility 
in terms of the right-hand-side of inequality Hll|) . This 
argument can be easily adopted for our system. If we 
assume that the main source of fluctuations of the in- 
stantaneous expression for the scattering function is the 
potential energy, we get 



6F,{k;t) 



dF,{k-t) 5V{Qi) 
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FIG. 5: Comparison of Xii^) (open cir- 

cles) with kBT^Xj^it) / cf^^ (closed circles) and 

{5Fs{k;t)&V{0)Y /{ksT^c^y') (open squares) at three 
representative temperatures: (a) T=1.00, (b) T=0.55 and 
(c) T=0.45. 



Eq. leads immediately to 



(14) 



dT 



(13) 



In Fig. we show the time dependence of both sides 
of the relation (|14|l at three representative temperatures. 
While the approximation 1)1411 is inaccurate, it becomes 
better, especially near the maximum value of X4(t), when 
the temperature approaches the mode coupling temper- 
ature. 
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FIG. 6: Temperature dependence of the maximum value of 
kBT^Xri^)/'^^ (closed circles) compared to that of the max- 
imum value of Xii't) (open circles). 



Figure El compares the temperature dependence of the 
maximum value of the right-hand-side of the relation H14(l 
to the temperature dependence of the maximum value of 
the four-point susceptibility. The difference between the 
two diminishes with decreasing temperature. However, it 
is clear that the temperature dependence of these quan- 
tities is different. In particular, one should be cautious 
when using the temperature dependence of the estimate 
(fnji to obtain the temperature dependence of the dy- 
namic correlation length, at least for temperatures above 
the mode coupling temperature. 

Finally, Fig. [7| compares the derivative of the loga- 
rithm of the four-point susceptibility, \nx4{t), with re- 
spect to the logarithm of time, Int and the derivative of 
the logarithm of the estimate lO, lnA;BT2xf.(t)/c^°* = 
2 In XT (^) + const., with respect to Int. We find that 
the estimate (|14|l approaches its maximum value with 
a power law exponent that is approximately equal to 1 
and thus somewhat higher than the power law exponent 
describing X4(t)'s approach to its maximum value. More- 
over, comparing Figs. [7^ and [3d we find that the differ- 
ence between exponents obtained from the estimate H14|l 
and from the four-point susceptibility does not seem to 
diminish with decreasing temperature. 



V. CONCLUSIONS 

We performed a quantitative analysis of the four-point 
susceptibility X4(t) of the Kob- Andersen Lennard- Jones 
binary mixture. We compared the results of Brownian 
dynamics computer simulations to predictions obtained 
from a recent re-formulation of the mode coupling theory. 
We did not compare computer simulation results to other 
approaches to glassy dynamics {e.g., to the predictions 
obtained using facilitated kinetic Ising models 0, IT^ ) 
because these other approaches are more concerned with 
a temperature range lower than the one accessible in our 
simulations. 

We found that some of the mode coupling predictions 
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FIG. 7: Comparison of the derivative din X4it)/ dint = 
[t/X4,{t))dx4,{t) /dt (open circles) and the derivative 
dlnkBT^XT{t)/c^y^ /dint = 2{t/xT{i))dxT{t)/dt (closed 
circles) at two temperatures: (a) T=0.55, and (b) T=0.45. 



agree with our simulation results. Most notably, the 
height of the maximum of the four-point susceptibility 
grows upon approaching the mode coupling temperature 
from above in the way predicted by the theory. More- 
over, the time at which x^if) reaches its maximum value 
has the same temperature dependence as the a relaxation 
time. For each of these two quantities, which are derived 
from a /owr-point function, the power law dependence on 
(T — Tc)/Tc is obeyed over the same temperature range 
as for quantities derived from two-point functions, e.g. 
the a relaxation time and the self-diffusion coefficient. 

On the other hand, we found some disagreement with 
the theory with regard to the power law dependence of 
Xi{t) on time. We found that upon approaching its max- 
imum value X4(t) grows with time with an effective ex- 
ponent of about 0.8. This exponent is quite a bit larger 
than the exponent predicted by the theory. It should be 
recalled that the exponent predicted by the theory is in 
turn larger than the one obtained by fitting a formula 
inspired by the mode coupling theory to the (two-point) 
self-intermediate scattering function. Hence, our under- 
standing of the connection between the time dependence 
of the four-point susceptibility and the time dependence 
of two-point functions seems incomplete. 

Finally, we tested the approximate estimate of the 
four-point susceptibility in terms of the temperature 
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derivative of the self-intermediate scattering function. 
We found that the estimate becomes accurate around 
the peak of Xiit) upon approaching the mode coupHng 
temperature. However, the temperature dependence of 
the maximum value of XiW is weaker than that of the 
approximate estimate. Moreover, the time dependence 
of X4:{t) differs from that of the estimate even near the 
mode coupling temperature. 

We should emphasize that all our computer simulation 
results were obtained using Brownian dynamics (and, 
therefore, NVT ensemble). It is possible that results ob- 



tained from Newtonian dynamics computer simulations 
would agree better (or worse) with theoretical predic- 
tions. 
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